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Summary 

Current  challenges  and  problems  pertaining  to  the  development  and  application  of  the 
DSMC  method  for  high-altitude  aerodynamics  are  discussed.  Attention  is  paid  to  issues 
related  to  the  efficiency  and  accuracy  of  the  method  in  the  near-continuum  regime,  as 
well  as  its  use  for  modeling  of  rarefied  flows  with  real  gas  effects.  Accurate  prediction 
of  force  and  heat  aerodynamic  characteristics  of  reentry  vehicles  and  spacecraft  requires 
comprehensive  investigation  of  hypersonic  flows  in  the  near-continuum  regime.  A  pow¬ 
erful  software  system  SMILE,  which  uses  the  direct  simulation  Monte  Carlo  method  as 
a  numerical  approach  and  is  built  basing  on  the  contemporary  knowledge  in  this  field, 
is  developed  to  solve  advanced  problems  of  high-altitude  aerothermodynamics.  Effective 
numerical  algorithms  and  physically  grounded  models  of  real  gas  effects  are  implemented 
in  SMILE.  Aerodynamics  of  promising  reentry  capsule  in  the  near-continuum  regime  is 
considered  as  an  example. 
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1  Introduction 


The  study  of  physical  phenomena  in  rarefied  nonequilibrium  flows  is  a  challenging 
problem  directly  related  to  the  development  of  new  aerospace  technologies.  Rarefied  gas 
dynamics,  that  deals  with  these  phenomena,  is  the  synthesis  of  a  great  number  of  funda¬ 
mental  problems  such  as  molecular  collision  dynamics  and  energy  transfer  phenomena  in 
collisions,  gas-surface  interactions,  condensation  and  evaporation,  plume  and  expansion 
flows,  and  many  others.  All  these  problems  are  in  close  connection  with  applied,  prac¬ 
tical  issues  that  can  be  conventionally  divided  into  two  groups.  The  first  group  covers 
the  questions  related  to  aerodynamic  calculations  of  hypersonic  flight  of  vehicles  at  high 
altitudes,  the  second  group  being  mainly  represented  by  the  problems  that  involve  the 
calculation  of  nozzle  flows  in  thrusters  and  jets  exhausting  into  vacuum  and  interacting 
with  the  surface  of  space  objects. 

Substantial  difficulties  arising  in  the  study  of  such  flows  are  caused  by  both  the  prob¬ 
lems  related  to  rarefaction  and  physico-chemical  effects.  It  is  commonly  known  that 
experimental  simulation  of  nonequilibrium  low-density  flows  is  rather  problematic  and 
expensive.  The  difficulties  of  experimental  modeling  have  stimulated  an  intense  develop¬ 
ment  of  various  approaches  for  numerical  simulation  of  these  flows.  Presently  there  are 
numerous  numerical  approaches  for  solving  the  problems  of  rarefied  gas  dynamics,  and 
the  choice  of  this  or  that  approach  depends  usually  on  the  flow  rarefaction,  the  problem 
dimension  and  the  presence  of  real  gas  effects. 

The  choice  of  the  numerical  approach  to  be  used  to  model  rarefied  nonequilibrium 
flows  greatly  depends  on  the  extent  of  flow  rarefaction.  For  near-continuum  flows,  it  is 
usually  sufficient  to  take  into  account  the  initial  effects  of  rarefaction  through  the  bound¬ 
ary  conditions  of  slip  velocity  and  temperature  jump  on  the  surface.  The  Navier-Stokes 
or  viscous  shock  layer  equations  are  commonly  used  with  these  boundary  conditions.  The 
Navier-Stokes  equations  can  be  derived  from  the  Boltzmann  equation  under  the  assump¬ 
tion  of  small  deviation  of  the  distribution  function  from  equilibrium.  Therefore,  they 
become  unsuitable  for  studying  rarefied  flows  where  the  distribution  function  becomes 
considerably  nonequilibrium. 

To  study  rarefied  flows  with  a  significant  degree  of  nonequilibrium,  the  Direct  Simu¬ 
lation  Monte  Carlo  (DSMC)  method  is  usually  employed.  This  method  has  become  the 
main  technique  for  studying  complex  multidimensional  rarefied  flows.  It  has  been  success¬ 
fully  applied  over  the  last  four  decades  to  model  various  flow  phenomena  and  gas  dynamic 
problems.  The  method  has  gradually  evolved  to  the  stage  where  its  application  to  cal¬ 
culate  complex  three-  dimensional  flows  is  almost  straightforward.  The  extended  area  of 
applicability  of  the  DSMC  method  is  from  the  near-continuum  regime  where  it  overlaps 
with  that  of  the  continuum  approaches,  to  the  free-molecular  regime.  In  practice,  the 
DSMC  method  is  computationally  intensive  compared  with  its  continuum  counterparts. 
However,  with  increasing  capabilities  of  modern  parallel  computers  this  method  acquires 
new  areas  of  application,  such  as  modeling  of  internal  and  external  near-continuum  flows, 
detailed  study  of  different  three-dimensional  problems  with  real  gas  effects,  and  others. 

The  Laboratory  of  Computational  Aerodynamics  from  the  Khristianovich  Institute  of 
Theoretical  and  Applied  Mechanics  of  the  Siberian  Branch  of  the  Russian  Academy  of 
Sciences  has  more  than  decade-long  experience  in  the  development  of  the  DSMC  method. 
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The  result  of  these  efforts  is  a  computational  system  called  SMILE  (Statistical  Modeling 
In  Low-Density  Environment)  capable  of  solving  a  very  wide  range  of  basic  and  applied 
problems.  The  SMILE  system  provides  a  complete  lifecycle  of  computations  starting  from 
a  geometry  model,  pre-processing,  going  through  the  computation  proper,  and  finishing 
with  post-processing  and  presentation  of  results.  All  SMILE  subsystems  have  a  Graphic 
User  Interface  (GUI),  which  makes  them  user-friendly  and  easy  to  use.  The  SMILE 
core  code  is  written  in  FORTRAN90  and  has  no  memory  limitations  specific  to  static 
FORTRAN  programs.  The  user  interface  of  the  system  is  written  in  C++  and  uses  a 
free  cross-platform  wxWidgets  GUI  library.  The  same  theoretical  basis  (majorant  colli¬ 
sion  frequency  schemes)  is  used  to  develop  more  advanced  software  systems  RGDAS  and 
SMILE++. 

The  principal  objectives  of  the  paper  are  to  analyze  the  modern  status  of  the  DSMC 
method,  outline  promising  directions  of  method  development,  discuss  problems  that  may 
be  encountered,  and  give  possible  ways  to  overcome  these  problems.  The  paper  discusses 
modeling  of  real  gas  effects,  the  numerical  accuracy  issues  for  the  DSMC  method,  and 
application  of  the  SMILE  software  system  for  studying  promising  capsule  aerodynamics. 


2  Conceptual  issues  of  the  DSMC  method 


The  DSMC  method  is  traditionally  considered  as  a  method  of  statistical  simulation 
of  the  behavior  of  a  great  number  of  simulated  gas  molecules.  Usually  the  number  of 
simulated  particles  is  large  enough  (~  105  —  108),  but  this  is  extremely  small  in  com¬ 
parison  with  the  number  of  real  molecules.  Each  simulated  particle  is  then  regarded  as 
representing  an  appropriate  number  of  actual  molecules,  Fnurn . 

The  main  principle  of  DSMC  is  the  splitting  of  continuous  process  of  molecular  motion 
and  collisions  into  two  successive  stages  at  the  time  step  At.  The  computational  domain 
is  divided  into  cells  of  size  Ax  such  that  the  variation  of  the  flow  parameters  in  every  cell 
is  small.  The  time  step  At  should  be  small  as  compared  with  the  mean  collision  time  T\. 
Free  motion  of  molecules  and  their  collisions  are  considered  successively  at  this  time  step 
At: 

1.  Collisions  of  particles  belonging  to  the  given  cell  in  each  cell  of  are  carried  out  inde¬ 
pendently  in  each  cell  of  physical  space,  i.e.  the  collisions  of  particles  in  the  neighboring 
cells  are  not  considered.  Since  the  distribution  function  variation  is  supposed  small  in  the 
cell,  when  a  pair  of  particles  for  collision  is  chosen  the  relative  distance  between  them  is 
not  taken  into  account.  The  post-collision  velocities  are  calculated  in  accordance  with  the 
conservation  laws  of  linear  momentum  and  energy. 

2.  All  molecules  located  in  the  computational  domain  are  displaced  by  the  distance 
determined  by  their  velocities  at  the  moment  and  by  the  time  step  At.  If  a  molecule  leaves 
the  computational  domain,  then  its  velocity  is  recomputed  according  to  the  boundary 
conditions.  At  the  same  step  At  new  particles  entering  the  computational  domain  are 
generated  in  accordance  with  the  distribution  function  specified  at  the  domain  boundaries. 

Thus,  the  following  principals  steps  are  specified  in  the  DSMC  procedure: 

•  entering  new  molecules 
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•  molecular  motion 

•  gas/surface  collisions 

•  indexing  of  molecules  over  collision  cells 

•  collisions 

•  sampling  of  macroparameters 

After  finishing  computation,  the  results  of  computations  are  processed  to  obtain  the 
flowfields  of  gas  dynamics  parameters,  total  and  distributed  surface  quantities.  The  main 
steps  of  the  DSMC  technique  are  discussed  below. 

The  state  of  each  simulated  particle  is  characterized  by  its  coordinate  f  and  velocity 
v.  The  state  of  the  whole  system  of  N  particles  is  described  by  a  6iV-dimensional  vector 
{R,  V}  =  {fi,  Vi, fjv,  hjv}.  The  evolution  of  such  a  system  can  be  represented  as  a 
jump-like  motion  of  a  point  in  the  6 AMlirriensional  phase  space.  The  DSMC  method  can 
be  then  treated  as  statistical  simulation  of  the  6-ZV-dimensional  random  jump- like  process. 

In  the  traditional  approach  Bird  [1]  of  constructing  numerical  schemes  of  the  DSMC 
method,  the  description  of  procedures  for  trajectory  simulation  of  the  random  process  is 
based  on  physical  concepts  of  rarefied  gas  and  on  physical  assumptions  that  create  the 
basis  for  the  phenomenological  derivation  of  the  Boltzmann  equation. 

A  different  approach  was  proposed  by  Ivanov  and  Rogasinsky  [2]  who  constructed 
numerical  schemes  for  simulating  a  6-ZV-dimensional,  continuous  in  time,  random  process 
of  spatially  nonuniform  evolution  of  a  system  of  N  particles.  Majorant  collision  frequency 
(MCF)  schemes  of  the  DSMC  method  were  derived  from  the  Kac  [3]  and  Leontovich 
[4]  master  kinetic  equations  (MKE).  This  is  a  linear  integro-differential  equation  that 
describes  the  time  behavior  of  an  iV-particle  gas  model  with  binary  collisions  on  the  level 
of  the  iV-particle  distribution  function  /at.  The  linear  MKE  may  be  transformed  into 
the  nonlinear  Boltzmann  equation  when  N  — »  oo  and  the  molecular  chaos  condition  is 
satisfied  (see,  e.g.,  [5]).  Since  a  finite  number  of  simulated  particles  is  used  in  the  DSMC 
simulations,  it  is  natural  to  directly  use  MKE  for  constructing  numerical  schemes  of  the 
DSMC  method.  The  procedure  of  constructing  MCF  schemes  is  described  in  detail  in 
Appendix  1. 

The  extension  of  the  majorant  collision  frequency  schemes  obtained  from  MKE  to  the 
case  of  multispecies  chemically  reacting  mixtures  is  straightforward.  It  implies  the  change 
in  cross-sections  of  the  corresponding  inelastic  processes  and,  hence,  the  change  in  the 
collision  algorithm  only  in  the  part  that  is  related  to  the  collision  mechanics. 

The  procedures  of  modeling  gas-surface  interaction  for  two  reflection  models  (specular- 
diffuse  and  Nocilla  models)  can  be  found  in  Appendix  2. 

Currently  used  numerical  schemes  of  the  DSMC  method  (NTC  [1],  NCT  [6],  MCF  [2, 
7])  are  very  close  both  in  terms  of  efficiency  and  numerical  implementation.  At  present, 
there  are  no  essential  internal  resources  for  increasing  the  efficiency  of  the  collision  stage 
of  the  DSMC  method.  The  main  effort  in  improving  the  method  efficiency  is  directed  to 
the  study  of  influence  of  the  grid  type,  the  use  of  multi-zone  approaches,  variable  time 
step,  grid  adaption  procedures,  etc  [8,9,38]. 
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3  Real  gas  effect  models  for  DSMC 


An  accurate  prediction  of  high-temperature  rarefied  flows,  such  as  those  behind  the 
shock  wave  formed  about  a  space  vehicle  at  high  altitudes,  requires  the  use  of  adequate 
models  of  physical  and  chemical  processes  -  so-called  real  gas  effects,  and  effective  numer¬ 
ical  procedures.  An  example  of  the  impact  of  these  processes  on  the  flow  is  given  in  Fig. 
1  where  the  translational  temperature  fields  about  the  Soyuz  capsule  at  an  altitude  of  85 
km  are  shown  for  nonreacting  and  reacting  gases  [9] . 

Two  most  important  effects  of  chemical  reactions  are  the  decrease  in  temperature  in 
the  shock  front  and  a  smaller  shock  stand-off  distance.  In  conventional  continuum  gas 
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Figure  1:  Translational  temperature  fields  about  a  reentry  capsule,  al¬ 
titude  85  km,  velocity  8  krn/s.  Nonreacting  (left)  and  reacting  (right) 

air. 


dynamics  the  real  gas  effects  are  usually  understood  as  such  high-temperature  phenomena 
as  molecular  vibration,  dissociation,  ionization,  surface  chemical  reactions,  and  radiation. 
In  describing  the  problems  of  rarefied  gas  dynamics,  with  a  typical  large,  often  of  the 
order  of  the  reference  flow  scale,  shock  wave  width  and  rarefaction  effects  exerting  the 
determining  effect  on  the  flow  structure,  it  is  convenient  to  consider  the  “real  gas  effects” 
in  a  wider  sense.  As  applied  to  kinetic  methods  based  on  the  microscopic  approach  and 
used  for  studying  RGD  problems,  it  seems  natural  to  relate  the  real  gas  effects  with  all 
phenomena  associated  with  molecular  collisions,  namely,  molecular  interaction  potential, 
rotational  degrees  of  freedom,  vibrational  degrees  of  freedom,  chemical  reactions  in  gas, 
ionization,  and  radiation. 

One  of  the  most  challenging  problems  remaining  in  terms  of  the  method  development 
and  improvement  is  related  to  the  need  to  effectively  and  reliably  simulate  processes  of 
energy  transfer  between  internal  and  translational  modes,  chemical  reactions,  ionization, 
and  radiation.  The  presence  of  one  or  more  of  these  processes  drastically  changes  flow 
properties  such  as  density  and  temperature,  and  the  development  of  adequate  models  is 
therefore  very  important.  A  number  of  DSMC  techniques  and  models  were  suggested 
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and  used  in  the  last  decade  that  cover  different  aspects  of  the  problem  (see,  for  example, 
Refs.  [10]- [14]  and  references  therein).  Currently,  there  is  a  lack  of  models  that  are 
general  enough  to  treat  various  molecular  interactions  and  processes,  sufficiently  accurate 
to  capture  complex  flow  physics,  easy  to  implement,  and  computationally  efficient  to  be 
applied  to  calculate  near-continuum  flows. 

In  what  follows  we  would  like  to  mention  the  collision  models  that  are  both  compu¬ 
tationally  efficient,  easy  to  implement,  and  sufficiently  general  to  be  applied  for  any  type 
of  interactions  between  atoms  and  diatomic  and  polyatomic  molecules.  The  Variable  Soft 
Sphere  (VSS)  model  [15]  is  one  of  the  most  widely  used  models  of  intermolecular  inter¬ 
action.  In  this  model  the  total  collision  cross-section  depends  on  the  relative  collision 
velocity,  and  two  parameters  of  the  model  are  determined  from  the  condition  of  coinci¬ 
dence  of  diffusion  and  viscous  cross-sections  of  the  VSS  and  the  inverse-power  potential 
models.  Even  though  the  model  does  not  include  the  attractive  part  of  the  potential,  it  is 
applicable  for  most  conditions  where  the  DSMC  is  used.  Rotational  energy  distribution 
is  conventionally  assumed  to  be  continuous,  which  makes  the  energy  transfer  algorithms 
simpler,  especially  for  polyatomic  molecules. 

The  discrete  description  was  suggested  in  [16]  for  different  types  of  polyatomic  molecules. 
The  traditional  Larsen-Borgnakke  (LB)  model  [17]  is  most  widely  used  in  the  DSMC 
method  to  describe  the  energy  exchange  between  the  translational  and  rotational  modes 
of  colliding  particles.  In  this  model,  the  energy  spectrum  of  rotational  mode  is  assumed 
continuous,  and  some  fraction  tp  =  1  /ZTi  of  the  total  number  of  collisions  takes  place 
with  the  energy  exchange  between  translational  and  rotational  degrees  of  freedom.  The 
post-collision  rotational  and  relative  translational  energies  are  simulated  in  accordance 
with  the  local  equilibrium  distribution  functions  and  are  proportional  to  the  number  of 
degrees  of  freedom  of  the  mode.  Several  improvements  were  suggested  for  the  LB  model 
to  make  the  rotational  relaxation  rate  in  the  DSMC  modeling  correspond  to  the  contin¬ 
uum  and  experimental  rates.  First,  temperature-dependent  rotational  collision  number 
was  suggested  [12],  Then,  a  correction  factor  was  derived  that  establishes  a  relation  of 
Zr  employed  in  the  DSMC  method  with  a  continuum  analog  used  in  the  Jeans  relax¬ 
ation  equation  [18].  Finally,  a  particle  selection  methodology  was  proposed  [19],  which 
prohibits  multiple  relaxation  events  during  a  single  collision,  and  matches  Jeans  equation 
for  general  gas  mixtures.  In  contrast  to  the  traditional  LB  model,  the  particle  selection 
methodology  prohibiting  multiple  relaxation  events  does  not  include  directly  rotational- 
rotational  energy  transfer.  Note  that  this  type  of  energy  transfer  has  not  been  studied 
separately  in  the  DSMC  method.  Since  the  translational-rotational  energy  transfer  is  fast, 
and  rotational  distribution  is  usually  less  important  for  chemical  reactions  and  radiative 
processes  than  vibrational  distribution,  accurate  modeling  of  the  translational-rotational 
energy  transfer  is  sufficient  for  most  cases. 

The  spacing  between  rotational  energy  levels  is  small,  and  the  continuum  model  of  ro¬ 
tational  mode  approximates  well  the  discrete  distribution  [20].  It  is  therefore  reasonable  to 
use  the  continuum  Larsen-Borgnakke  model  with  the  temperature-dependent  relaxation 
number,  the  correction  factor  [18]  and  a  selection  methodology  that  prohibits  multiple 
relaxation  events.  Two  models  for  the  vibrational  energy  distributions,  continuous  and 
discrete,  are  traditionally  used  in  the  DSMC  method.  Generally,  the  continuity  of  the 
vibrational  energy  mode  may  be  a  too  rough  allowance,  since  the  vibrational  spectrum  of 
real  molecules  is  characterized  by  large  gaps  between  the  neighboring  energy  levels.  Be- 
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sides  the  physical  reason,  there  are  several  numerical  problems  connected  with  the  use  of 
the  continuum  vibrational  energy  spectra.  These  are  the  need  for  special  cut-off  parame¬ 
ters  for  the  LB  model  when  the  number  of  vibrational  degrees  of  freedom  is  less  than  two, 
additional  assumptions  that  have  to  be  used  in  the  LB  collision  energy  transfer  modeling, 
and  an  inaccurate  shape  of  the  internal  energy  distribution  functions  at  equilibrium  [20]. 
The  last  point  is  important  for  reactions  whose  cross-sections  directly  depend  on  vibra¬ 
tional  states  (such  as  vibrationally  favored  reactions).  Also,  internal  energy  distribution 
strongly  affects  the  flow  radiation.  The  discrete  internal  energy  model  does  not  have  these 
drawbacks,  but  requires  a  special  correction  procedure  to  be  used  for  the  chemical  reac¬ 
tion  rates  to  follow  the  experimental  rates  expressed  in  the  Arrhenius  form.  Many  discrete 
models  for  the  translation- vibration  energy  transfer  were  suggested  in  the  literature.  They 
differ  considerably  in  the  way  of  determining  the  vibrationally  inelastic  cross-sections.  A 
discrete  version  of  the  LB  model  is  presented  in  [21].  In  the  model  [22],  the  probability 
of  VT  transitions  is  determined  using  the  inverse  Laplace  transformation  for  a  simple 
harmonic  oscillator  from  the  temperature-dependent  relaxation  rate  [23].  The  approach 
[24]  is  based  on  the  use  of  the  information  theory  for  determining  post-collision  states.  In 
the  model  [25]  the  cross-sections  for  vibration-translation  transitions  are  used,  obtained 
for  an  anharmonic  oscillator  from  a  quasi- classical  approximation  of  the  scattering  theory. 

The  vibration- vibration  energy  transfer  was  also  modeled  in  the  DSMC  method.  We 
mention  here  the  model  [26]  where  both  vibration-translation  and  vibration-vibration 
transitions  for  a  simple  harmonic  oscillator  have  been  considered.  The  vibration- vibration 
energy  exchange  model  based  on  the  quasiclassical  approach  was  developed  in  [25].  This 
energy  transfer  process  was  found  to  be  important  for  the  vibrational  populations  in 
hypersonic  flows  at  high  altitudes.  All  these  and  many  other  models  have  been  developed 
only  for  diatomic  molecules.  Because  of  this,  it  seems  most  reasonable  at  present  to 
use  the  LB  model  with  particle  selection.  Temperature-dependent  vibrational  relaxation 
rate  has  to  be  specified  using  vibrational  mode  characteristic  temperatures  of  polyatomic 
molecules.  A  correction  factor  [27]  has  to  be  used  that  enables  one  to  match  the  vibrational 
relaxation  rates  in  the  DSMC  modeling  to  the  relaxation  rate  given  by  the  Landau-Teller 
equation.  Note  that  the  LB  model  was  not  yet  applied  to  vibration-vibration  energy 
transfer.  Such  an  application  would  require  to  specify  the  corresponding  collision  numbers 
that  may  not  be  known  for  many  molecular  systems. 

The  DSMC  method  has  been  used  for  calculating  rarefied  hypersonic  chemically  re¬ 
acting  flows  for  more  than  two  decades.  The  major  problem  in  generating  the  model 
of  chemical  reactions  is  the  determination  of  energy-dependent  cross-sections  of  chemi¬ 
cal  reactions.  Note  here  that  the  temperature  dependences  of  chemical  reaction  rates, 
traditional  for  the  continuum  approach,  cannot  be  applied  to  the  DSMC  method,  and 
energy-dependent  reaction  cross-sections  have  to  be  used.  As  the  exact  expressions  for 
these  cross-sections  are  not  available  now  even  for  simple  dissociation  reactions  in  air, 
some  approximations  have  to  be  employed.  The  first  simplistic  models  have  been  replaced 
by  the  total  collision  energy  (TCE)  model  [28]  built  on  the  basis  of  collision  theory  for 
chemical  reactions.  This  model  is  efficient  and  is  used  presently  for  calculating  two-  and 
three-dimensional  flows.  It  employs  the  major  assumption  that  the  reaction  probability 
depends  on  the  total  collision  energy.  A  specific  form  of  this  dependence  is  assumed  that 
allows  analytic  determination  of  unknown  coefficients.  The  derivation  uses  the  reaction 
rate  coefficient  in  Arrhenius  form  and  the  Hinshelwood  expression  for  the  equilibrium  en- 
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ergy  distribution  function  that  implies  continuous  distribution  of  internal  modes  [28].  The 
problems  inherent  in  continuous  vibrational  energy  models  that  were  discussed  above  lead 
to  the  conclusion  that  the  use  of  discrete  models  in  the  DSMC  method  is  preferable.  One 
problem  related  to  the  use  of  discrete  energy  levels  is  that  the  existing  chemical  reaction 
models  derived  for  continuous  internal  energies  need  special  modifications  to  be  applied 
for  discrete  energies.  In  this  case,  a  correction  procedure  is  necessary  for  the  reaction  rates 
in  the  DSMC  method  match  available  experimental  data.  Such  a  correction  procedure 
was  developed  in  [20].  It  utilizes  a  Monte  Carlo  approach  to  find  modified  constants  in 
the  reaction  cross-section  energy  dependence,  the  use  of  which  enables  one  to  match  the 
original  reaction  rates  at  equilibrium. 

The  procedure  is  applied  to  the  total  collision  energy  model,  but  is  generally  applicable 
to  any  reaction  model  based  on  the  collision  theory  for  chemical  reactions.  Several  models 
that  include  a  direct  dependence  of  the  reaction  cross-sections  on  the  collider  vibrational 
energy  (so-called  vibration-dissociation  coupling)  were  proposed  earlier  [10,  13,  14],  The 
deficiency  of  these  models  is  that  they  imply  some  fixed  or  variable  parameters  that 
determine  the  extent  of  vibrational  favoring.  These  parameters  are  extracted  or  verified 
through  comparison  with  available  experimental  data.  Note  that  some  reactions  may  not 
have  vibrational  favoring,  and  therefore  the  model  [28]  may  turn  to  be  reasonable  [29]. 


4  DSMC  versus  continuum  CFD 


The  areas  where  the  continuum  and  kinetic  approaches  are  currently  applied  over¬ 
lap  in  the  regime  of  low  Knudsen  numbers.  The  recent  experimental  and  computational 
activity  related  to  accurate  prediction  of  laminar  flow  separation  in  the  near-continuum 
regime  shows  that  even  for  non-reacting  gases  there  are  still  problems  of  validation  and 
verification  of  numerical  approaches  that  need  to  be  resolved.  Numerical  analysis  of  such 
flows  is  traditionally  performed  using  Navier-Stokes  (NS)  equations,  with  the  initial  ef¬ 
fects  of  rarefaction  taken  into  account  through  the  slip  velocity  and  temperature  jump. 
The  use  of  slip  conditions  is  caused  by  the  fact  that  rarefaction  effects  can  be  observed 
for  slender  bodies  even  for  rather  high  Reynolds  numbers  (Re  =  20,000  —  30,000).  The 
state-of-the-art  algorithms  of  the  DSMC  method  (adaptive  grids,  variable  time  step,  etc.) 
allow  one  to  simulate  flows  at  such  high  Reynolds  numbers.  This  makes  realistic  the 
study  of  the  applicability  area  of  the  continuum  approach  to  model  laminar  separated 
flows.  The  study  was  started  in  [30]  and  continued  in  [31]  where  a  detailed  numerical 
modeling  of  an  axisymmetric  shock  wave  /  laminar  boundary  layer  interaction  was  per¬ 
formed  by  the  continuum  (NS)  and  kinetic  (DSMC)  approaches  for  the  ONER  A  R5Ch 
wind  tunnel  conditions.  The  DSMC  results  showed  the  presence  of  a  high  slip  velocity 
and  temperature  jump  near  the  wall.  The  application  of  slip  conditions  for  the  NS  solver 
has  a  substantial  effect  on  the  entire  flow  held  and  allows  one  to  decrease  the  difference 
between  the  continuum  and  DSMC  results.  The  left  part  of  Fig.  2  shows  that  the  profile 
of  the  gas  velocity  at  the  wall  Ug  obtained  by  the  DSMC  method  is  in  excellent  agreement 
with  the  Ug  profile  calculated  from  the  NS  slip  velocity  using  the  solution  for  the  velocity 
distribution  across  the  Knudsen  layer  [32],  However,  the  NS  solver  with  slip  conditions 
predicts  a  larger  length  of  the  separation  region  than  the  DSMC  method.  To  eliminate 
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U/U«  Y/L  Y/L 


Figure  2:  Slip  and  gas  NS  velocities  and  gas  DSMC  velocity  at  the  wall  (left).  Density 
(center)  and  temperature  (right)  profiles  for  a  monatomic  gas. 


the  effect  caused  by  different  descriptions  of  translational/rotational  energy  exchange  for 
both  methods,  NS  and  DSMC  simulations  of  a  laminar  separated  flow  of  a  monatomic  gas 
(argon,  M  =  10)  were  performed  in  [32],  The  central  part  of  Fig.  2  shows  the  comparison 
of  density  profiles  in  the  cross-section  X/L  =  0.3.  The  use  of  slip  conditions  for  the  NS 
solver  reduces  the  slope  of  the  leading-edge  shock  wave  and  generally  gives  a  good  agree¬ 
ment  with  the  DSMC  results.  The  comparison  of  temperature  profiles  (right  part  of  Fig. 
2)  shows  that  the  difference  between  the  no-slip  continuum  and  kinetic  results  reaches 
40%.  The  use  of  the  slip  conditions  allows  one  to  obtain  the  temperature  near  the  body 
surface  that  is  very  close  to  the  DSMC  values,  and  also  decrease  the  difference  inside  the 
boundary  layer.  Still,  the  temperatures  inside  the  boundary  layer  for  DSMC  and  NS  with 
slip  conditions  differ  by  about  20%.  Figure  3  shows  the  distributed  aerothermodynamic 


X/L  X/L  X/L 


Figure  3:  Distributed  surface  characteristics  for  monatomic  gas.  Pressure  coefficient  (left), 
skin  friction  coefficient  (center),  and  Stanton  number  (right) 


characteristics  obtained  using  the  continuum  and  kinetic  approaches.  The  differences  in 
the  pressure  and  skin  friction  coefficients  and  the  Stanton  number  are  mainly  due  to  the 
following  two  factors.  First,  the  length  of  the  separation  region  is  different.  An  earlier 
flow  separation  in  NS  results  leads  to  a  lower  pressure  peak  and  smaller  heat  flux  on  the 
flare  surface.  Second,  there  is  a  strong  flow  nonequilibrium  near  the  leading  edge.  This 
is  why  the  NS  solver  with  either  no-slip  or  slip  conditions  predicts  a  significantly  larger 
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heat  flux  (up  to  a  factor  of  four)  than  the  DSMC  method  in  the  vicinity  of  the  leading 
edge.  This  also  leads  to  a  significant  decrease  in  the  flow  temperature  downstream  and 
an  earlier  separation  of  the  boundary  layer.  The  studies  conducted  for  a  hypersonic  flow 
around  a  hollow-cylinder  flare  clearly  showed  that  the  Navier-Stokes  equations,  even  with 
slip  conditions,  are  inapplicable  in  the  vicinity  of  the  leading  edge.  The  flow  near  the 
leading  edge  affects  the  flow  structure  downstream,  for  example,  leads  to  an  increase  in 
the  separation  region  length. 


5  Numerical  accuracy  of  the  DSMC  method 


An  important  problem  that  arises  when  modeling  near-continuum  flows  is  the  eval¬ 
uation  of  the  numerical  accuracy.  Applying  the  DSMC  method  to  compute  transitional 
flows  (. Kn  ~  1),  one  can  usually  perform  additional  computations  on  a  finer  grid  and 
with  a  larger  number  of  simulated  molecules  in  order  to  analyze  the  numerical  accuracy. 
For  a  near-continuum  regime,  the  computations  are  often  conducted  using  all  available 
computer  resources,  and  no  additional  accuracy-establishing  computations  are  possible  in 
this  case.  The  question  is  therefore  what  is  the  accuracy  of  results  obtained,  and  how 
far  they  are  from  the  solution  of  the  Boltzmann  equation.  This  question  is  especially 
important  for  the  near-continuum  flows. 

The  sources  of  numerical  error  have  been  analyzed  in  the  continuum  CFD  for  many 
years.  There  are  many  techniques  available  such  as  the  analysis  of  truncation  errors, 
the  grid  convergence  study,  etc.,  that  enable  one  to  estimate  the  accuracy  of  the  results. 
There  is  nothing  similar  for  the  DSMC  method.  The  accuracy  analysis  in  DSMC  is 
complicated  by  the  presence  of  different  errors  that  may  be  classified  as  follows:  (1) 
statistical  errors;  (2)  spatial  resolution  and  time  step  errors;  (3)  errors  due  to  the  finite 
number  of  simulated  molecules.  The  complete  analysis  of  the  accuracy  should  include  the 
derivation  of  the  governing  equation  that  describes  the  modeling  process,  and  establishing 
of  the  connection  between  this  equation  and  the  Boltzmann  equation.  This  connection 
can  hardly  be  established  for  the  NTC  and  NCT  schemes  because  these  schemes  were 
formulated  at  the  phenomenological  level. 

The  MFS  scheme  was  derived  from  MKE,  which  describes  the  behavior  of  an  Y-particle 
gas  model  with  binary  collisions.  This  equation  may  be  transformed  to  the  Boltzmann 
equation  when  N  — »  oo  and  the  molecular  chaos  condition  is  satisfied.  Let  us  examine 
the  molecular  chaos  requirement  in  more  detail.  A  finite  number  of  molecules  N  is  always 
used  in  DSMC  modeling.  Therefore,  there  is  a  statistical  dependence  between  simulated 
particles  [33].  The  magnitude  of  this  dependence  depends  on  the  total  number  of  simulated 
molecules  in  the  system,  and  generally  decreases  when  the  number  of  molecules  increases. 
It  is  necessary  to  estimate  the  level  of  statistical  dependence  and  its  contribution  to  the 
results  of  DSMC  computation.  A  significant  level  of  statistical  dependence,  or  particle 
correlations,  means  that  the  molecular  chaos  hypothesis,  used  in  the  Boltzmann  equation, 
is  no  longer  valid.  The  molecular  chaos  hypothesis  may  be  written  in  terms  of  one-  and 
two-particle  distribution  functions,  f\  and  fi  ,  as 

f2(vA,VB)  =  fiMfiM, 
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where  vA  and  uB  are  the  velocities  of  the  two  collision  partners,  A  and  B.  For  some 
function  of  velocity,  h  =  h(u),  we  can  write 

/oo  /»oo  /»oo 

f2{vA,vB)h(vA)h(vB)dvAdvB=  /  fi{vA)h{yA)dvA  /  fi{yB)h(vB)dvB 

00  J—  oo  J  —  OO 

This  expression  may  be  rewritten  as 

(h(isA)h(vB))  =  (h{yA))(h{yB)) 

where  (...)  denotes  averaging.  As  was  mentioned  above,  the  molecular  chaos  hypothesis  is 
not  strictly  valid  for  a  finite  number  of  particles,  which  means  that  the  above  expression 
does  not  hold.  Let  us  introduce  the  following  correlator  G2, 

=  {h(vA)h(vB))  _  1 
2  (h(uA))(h(uB)) 

Obviously,  this  correlator  may  serve  as  an  indicator  of  the  presence  of  molecular  corre¬ 
lations.  The  smaller  the  correlations,  the  closer  G2  is  to  0.  The  value  of  G2  may  be 
directly  calculated  in  the  DSMC  method  over  all  cells  in  the  computational  domain,  and 
h(y)  =  V2  is  usually  used  for  this  purpose  [34],  [35].  An  additional  criterion  that  allows 
for  a  practical  verification  of  the  presence  of  the  statistical  dependence  between  simu¬ 
lated  particles  is  the  relative  number  of  repeated  collisions  [34],  Repeated  collisions  are 
collisions  between  the  same  pair  of  particles  during  their  lifetime  in  the  computational 
domain.  The  number  of  repeated  collisions  is  directly  related  to  the  number  of  particles 
N\  in  a  volume  with  the  linear  size  equal  to  the  local  mean  free  path.  If  N\  >  1,  one  can 
say  that  the  simulation  results  are  close  to  the  solution  of  the  Boltzmann  equation.  Let 


Figure  4:  ID  shock  wave,  monatomic  gas.  Density  n  and  parallel  temperature  Tx  (left). 
Correlator  G2( V2)  and  repeated  collision  fraction  (right). 


us  consider  a  1 D  shock  wave  structure  as  an  example  where  the  G2  correlator  and  the 
repeated  collision  number  are  used.  A  monatomic  gas  of  hard  sphere  molecules  is  used 
with  M  =  8.  The  computational  domain  was  15  upstream  mean  free  paths.  Three  cases 
were  considered  here.  20  particles  in  the  upstream  cells  was  used  in  Case  1  (baseline  case), 
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with  the  cell  size  L  of  0.25  of  the  upstream  mean  free  path  A.  For  Case  2,  2  particles  per 
upstream  cell  was  used  with  L  =  0.25A.  The  number  of  particles  for  case  3  was  the  same 
as  in  Case  1,  but  L  =  0.00625A.  The  average  number  of  particles  in  the  upstream  cells 
was  therefore  0.5  in  this  case. 


The  density  and  parallel  temperature  distributions  inside  the  shock  wave  are  shown 
in  Fig.  4  (left).  The  profiles  are  very  similar  for  Cases  1  and  3,  and  somewhat  different 
for  Case  2.  The  G2  correlator  and  repeated  collision  fraction  are  given  in  Fig.  4  (right). 
For  Case  1,  the  values  of  G2  are  close  to  zero,  being  less  than  0.01  even  in  the  region 
of  strong  nonequilibrium,  and  the  number  of  repeated  collisions  is  not  greater  than  2 
percent  of  the  total  number  of  collisions.  The  values  of  the  correlator  and  repeated 
collision  fraction  increase  about  ten  times  for  Case  2  compared  to  case  1,  even  though  the 
collision  frequency  is  the  same  since  the  MFS  scheme  is  used.  The  increase  in  the  particle 
statistical  correlations  is  the  reason  for  the  difference  in  the  macroparameters  observed 
in  Fig.  4  (left).  Obviously,  the  difference  is  due  to  the  decrease  in  the  total  number  of 
simulated  particles  in  the  computational  domain.  We  can  therefore  state  that  an  accurate 
calculation  of  the  collision  frequency  in  a  cell  is  not  sufficient  for  obtaining  an  accurate 
result.  An  interesting  fact  is  that  the  reduction  of  the  cell  size  while  preserving  the  total 
number  of  molecules  (Cases  1  and  3)  does  not  change  macroparameters  and  only  slightly 
changes  the  correlator  and  the  number  of  repeated  collisions.  This  shows  that  the  number 
of  particles  in  cell  is  not  a  determining  factor  for  MFS.  The  important  parameters  are  the 
total  number  of  particles  in  the  system  and  N\. 


This  fact  is  illustrated  below  for  a  flow  over  a  cylinder  in  CO2/N2  mixture  [36].  A 
free-stream  velocity  of  7400  m/s,  temperature  of  137.4  K  and  Knudsen  number  Kn  =  0.01 
were  taken.  The  mole  fractions  of  CO2  and  N2  were  assumed  to  be  0.9537  and  0.0463, 
respectively.  A  diffuse  reflection  with  complete  energy  accommodation  was  employed 
at  the  surface  with  temperature  value  of  1000  K.  The  computations  were  performed  for 
different  numbers  of  simulated  molecules  in  the  computational  domain  in  order  to  examine 
the  impact  of  N\  on  results.  An  adaptive  grid  refinement  technique  was  used  to  satisfy 
the  Knc  >  1  condition  in  the  entire  computational  domain.  The  profiles  of  translational 
and  vibrational  temperatures  along  the  stagnation  line  are  plotted  in  Fig.  5  (X  =  0 
corresponds  to  the  stagnation  point).  The  results  are  shown  for  different  values  of  N\ 
in  the  free  stream  which  corresponds  to  different  values  of  the  total  number  of  simulated 
particles.  The  N\  distribution  along  the  stagnation  line  is  also  presented  here.  The  local 
N\  is  not  less  than  1  for  the  free  stream  N\  =  20,  and  this  case  is  therefore  considered 
to  be  the  reference.  It  is  seen  that  for  N\  =  2.5  in  the  free  stream  the  results  agree  with 
the  reference  case  inside  the  shock  wave,  but  are  somewhat  different  near  the  stagnation 
point.  This  is  because  N\  in  this  region  decreases  sharply  and  becomes  less  than  1.  When 
the  free  stream  N\  is  further  decreased,  the  results  are  relatively  close  to  the  reference 
case  inside  the  shock,  but  significantly  differ  inside  the  boundary  layer.  It  should  be  noted 
that  the  impact  of  N\  is  more  pronounced  for  the  vibrational  temperature  than  for  the 
translational  one.  The  results  presented  show  that  the  condition  N\  >  1  must  be  satisfied 
in  DSMC  simulations  for  the  impact  of  particle  correlations  be  insignificant. 
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Figure  5:  Translational  Tt  (left)  and  vibrational  Tv  (center)  temperatures  and  N\  (right) 
along  the  stagnation  line  for  different  free  stream  N\  values. 


6  High-altitude  aerothermodynamics  of  a  promising 
capsule 


PPTS  (Prospective  Piloted  Transport  System),  unofficially  called  Rus,  is  a  project  be¬ 
ing  undertaken  by  the  Russian  Federal  Space  Agency  to  develop  a  new-generation  manned 
spacecraft.  Its  official  name  is  Pilotiruemyi  Transportny  Korabl  Novogo  Pokoleniya  or 
PTK  NP  meaning  New  Generation  Piloted  Transport  Ship.  The  goal  of  the  project  is  to 
develop  a  new-generation  spacecraft  to  replace  the  aging  Soynz  which  was  developed  by 
the  former  Soviet  Union  more  than  forty  years  ago. 

At  the  initial  stage  of  spacecraft  design,  it  is  necessary  to  study  its  aerothcrmodynamic 
characteristics  in  a  wide  range  of  free-stream  parameters.  During  de-orbiting,  the  vehicle 
passes  through  the  free-molecular  flow,  then  through  the  transitional  zone,  and  the  flight 
is  finalized  in  the  continuum  flow.  The  loads  on  the  vehicle  during  its  descent  are  consid¬ 
erably  changing.  The  efficiency  of  control  surfaces  is  also  changing.  To  predict  the  zone 
of  spacecraft  landing,  one  should  know  the  aerodynamic  characteristics  of  the  vehicle  at 
all  stages  of  its  descent.  The  results  of  statistical  simulation  of  aerothermodynamics  of  a 
promising  space  capsule  at  altitudes  from  120  to  60  km  are  presented  below  (for  details, 
see  [47]). 

The  aerothermodynamic  parameters  of  the  space  vehicle  at  altitudes  from  120  to  60 
km  were  studied  by  the  local  bridging  method.  Standard  atmosphere  parameters  were 
taken  for  computations  from  GOST  4401-81.  The  DSMC  method  was  used  to  model  the 
flow  at  several  altitudes  with  the  parameters  listed  in  Table  1. 


H,  km 

75 

80 

85 

90 

100 

110 

Moo 

25.9 

26.5 

27.2 

27.3 

26.3 

22.5 

Too,  K 

208.4 

198.6 

189.0 

186.7 

196.6 

255.5 

Poo,  kg/m 3  -106 

39.9 

18.5 

8.21 

3.42 

0.56 

0.093 

Klloo 

0.38  •  HT3 

0.8  •  HT3 

0.18  •  HT2 

0.43  •  10~2 

0.27  •  10_1 

0.17 

Table  1:  Parameters  of  the  atmosphere  (GOST  4401-81) 


The  space  vehicle  considered  is  shaped  similar  to  the  Apollo  capsule.  The  nose  part  is 
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spherical,  and  the  rear  part  is  beveled  at  an  angle  of  20°.  The  vehicle  geometry  is  shown  in 
Fig.  6.  A  qualitative  difference  from  the  Apollo  capsule  is  the  presence  of  trimming  flaps. 
One  of  the  objectives  of  the  present  activities  was  to  determine  the  changes  in  efficiency 
of  these  flaps  with  variations  of  the  angle  of  attack  and  with  allowance  for  real  gas  effects. 
The  temperature  of  the  entire  body  was  assumed  to  be  constant  and  equal  to  1000  K. 

The  influence  of  chemical  reactions  proceeding  in  the  flow  on  the  efficiency  of  control 
surfaces  was  determined  in  an  axisymmetric  formulation.  As  the  axisymmetric  formu¬ 
lation  allows  us  to  obtain  the  aerodynamic  characteristics  of  only  those  elements  that 
are  obtained  by  means  of  rotation  by  360°,  but  the  trimming  flaps  have  finite  size,  their 
contribution  to  the  aerodynamic  parameters  was  taken  into  account  with  the  following 
procedure:  the  aerodynamic  loads  (Cx  and  Cy)  on  a  circular  flap  were  identified;  then 
the  fraction  of  the  loads  on  the  flap  in  accordance  with  its  angular  size  was  calculated; 
finally,  the  pitching  moment  generated  by  the  flaps  was  found  from  the  loads  obtained 
(Cx  flap  and  Cy  flap)- 

The  DSMC  computations  were  performed  with  the  use  of  the  SMILE  software  system 
[7] .  The  internal  energy  exchange  was  modeled  in  accordance  with  the  Larsen-Borgnakke 
model.  Chemical  reactions  that  occur  during  particle  collisions  were  also  taken  into  ac¬ 
count.  Diffuse  reflection  from  the  surface  with  complete  accommodation  of  energy  was 
assumed. 

The  computations  were  performed  on  clusters  of  the  Interdepartmental  Supercomputer 
Center  (Moscow,  Russia)  and  of  the  Khristianovich  Institute  of  Theoretical  and  Applied 
Mechanics  of  the  Siberian  Branch  of  the  Russian  Academy  of  Sciences  (Novosibirsk,  Rus¬ 
sia).  Up  to  128  processors  were  used.  The  computations  at  an  altitude  of  75  km  required 
approximately  7000  processor-hours. 

The  local  bridging  method  used  in  this  work  is  based  on  the  technique  proposed  in  [46] . 
The  influence  of  the  flow  on  an  elementary  area  is  described  by  the  following  formulas: 

P  =  P0  +  PlWn  +  P2W2n  T  =  T0Wt  +  TiWnWt  (1) 

Here  P  and  r  are  the  pressure  and  friction  coefficients  of  aerodynamic  forces  normalized 
to  the  dynamic  pressure  q  =  ,  w  is  the  normalized  free-stream  velocity  vector, 

and  wn  and  wt  are  the  normal  and  tangential  components  of  the  velocity  vector  to  the 
surface.  The  parameters  P0,  Pi,  P2,  r0,  and  iq  are  calculated  using  the  following  formulas: 


P0  =  P«  +  (Pt  -  P^Fpo  Pi  =  P(mFp2  P2  =  P\d  +  (. Plm  -  P2d)FP2  (2) 

t0  =  rlmFr  o  Ti  =  r{mFTi  (3) 

(n  is  the  normal  vector  to  the  surface  and  t  is  the  tangential  vector).  These  vectors  lie  in 
one  plane. 

The  bridging  coefficients  FP0,  FP1,  FP2,  Frl,  and  Frl  are  used  to  take  into  account 
the  influence  of  the  free-molecular  or  continuum  flow  regime  under  particular  conditions. 

The  formulas  for  calculating  the  bridging  functions  FP  and  Fr  were  derived  serni- 
empirically  in  [46]  Here  we  give  their  final  form: 
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Figure  6:  Space  vehicle  geometry 


a\/Reo  +  exp(— 5Re0)  ’ 
(7-1)\/G;  +  M-1v/2(7-1) 
(0.56 +  1.2t^)(M  + 2.15) 


where 


b  =  0.35  +  0.005M 


(4) 

(5) 


Fpi  =  FP2  =  exp  (-(0.125  +  0.078^)Re0  •  (6) 


a  i 


Fto  =  [aiRe0  +  exp(-6iRe0)] 


7  ~  1 
2 


M(0.208  + 0.341  tw) 


■3/4 

5 

-3/4 


where 

&i  =  0.213-0.133^ 


(7) 

(8) 


Fti  =  [0.145/2  +  exp  (7.2  •  10_3i?  —  1.6  •  10-5i?2)]  1/2  ,  where  (9) 

R  =  (0.75 tw  +  0.25)_2/3  Re0  •  10"2-4(1-sinai)3  (10) 

Figure  7  shows  the  pressure  fields  calculated  for  non-deflected  control  flaps  and  for  flaps 
deflected  by  30°.  The  calculations  were  performed  for  a  gas  without  chemical  reactions.  As 
it  could  be  expected,  the  flap  influence  extends  to  a  fairly  small  distance  upstream,  but  the 
field  in  the  wake  behind  the  vehicle  is  noticeably  changed.  Thus,  for  instance,  the  pressure 
in  the  middle  of  the  vehicle  base  area  is  approximately  90  Pa  for  the  non-deflected  flap 
and  about  50  Pa  for  the  flap  deflected  by  30°.  In  a  hypersonic  flow  (M  =  25.9),  the  base 
pressure  does  not  exert  any  significant  influence  on  the  values  of  the  total  aerodynamic 
coefficients.  Flap  deflection  increases  the  drag  coefficient  Cp>  from  1.50  to  1.54.  The  lift 
coefficient  is  Cp  =  0.103  for  the  flap  deflected  by  30°. 

The  effect  of  chemical  reactions  on  the  flow  field  is  illustrated  in  Fig.  8,  which  shows 
the  pressure  fields  around  the  space  vehicle  obtained  in  a  chemically  reacting  flow  (upper 
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P,  Pa 
3 . 254e+00 
6 . 779e+00 
1 . 061e+01 
1 . 479e+01 
1. 937e+01 
2 . 440e+01 
2 . 997e+01 
3. 616e+01 
4 . 307e+01 
5. 084e+01 
5. 966e+01 
6. 973e+01 
8 . 135e+01 
9 . 491e+01 
1 . 109e+02 
1 . 302e+02 
1 . 537e+02 
1 . 830e+02 
2 . 208e+02 
2 . 712e+02 
3 . 417e+02 
4 . 474e+02 
6 . 237e+02 
9 . 762e+02 


Figure  7:  Pressure  fields.  Altitude  75  km.  Chemically  inert  gas.  Non-deflected  flap 
(upper)  and  flap  deflected  by  30° (downer). 


figure)  and  in  a  chemically  inert  gas  (lower  figure)  for  the  flap  deflected  by  30°.  The 
calculations  were  performed  for  an  altitude  of  75  km.  It  is  seen  that  the  flow  fields  are 
drastically  different  in  the  entire  computational  domain.  Indeed,  as  the  major  part  of  the 
free-stream  energy  is  spent  on  dissociation  of  diatomic  molecules,  the  flow  temperature  in 
the  case  of  a  chemically  reacting  gas  is  substantially  lower  than  the  temperature  obtained 
with  the  chemical  reactions  being  ignored.  In  addition  to  temperature  fields,  the  Mach 
number  fields  are  also  essentially  different.  Therefore,  we  can  conclude  that  the  entire 
flow  structure  and  the  shock  waves  near  the  vehicle  are  formed  in  a  different  manner. 

p,  Pa 

I  3.246e+00 
U  6.762e+00 
1 .058e+01 
1.475O+01 
1.932O+01 
2.434e+01 
2.989e+01 
3.606e+01 
4.296e+01 
5.071  e+01 
5.950e+01 
6.955e+01 
8.114e+01 
9.467e+01 
1 1 .1 06e+02 
1 1 ,298e+02 

1 1.533O+02 

1 1 ,826e+02 
I  2.202e+02 
H  2.705e+02 
3.408e+02 
4.463e+02 
6.221  e+02 
9.737e+02 


Figure  8:  Pressure  fields.  Altitude  75  km.  Flap  deflection  angle  30°.  Chemically  reacting 
gas  (upper  figure)  and  chemically  inert  gas  (lower  figure). 


Figure  9  shows  the  distributions  of  temperature  and  pressure  normalized  to  the  free- 
stream  pressure  over  the  stagnation  lines  at  altitudes  of  75,  85,  and  90  km.  The  plots 
show  the  stagnation  lines  obtained  for  chemically  reacting  and  chemically  inert  gases. 
The  plots  on  the  left  show  the  free-stream  parameters,  and  the  plots  on  the  right  show 
the  situation  with  the  body  located  at  X  =  0.  It  is  seen  that  the  temperature  of  the 
flow  behind  the  bow  shock  wave  is  substantially  lower  if  chemical  reactions  are  taken 
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into  account.  In  addition,  the  stand-off  distance  of  the  bow  shock  wave  is  changed.  The 
influence  of  chemical  reactions  becomes  attenuated  with  increasing  flight  altitude.  Thus, 
for  instance,  the  difference  in  the  shock-wave  stand-off  distance  is  0.187  m  at  an  altitude 
of  75  km,  0.133  m  at  85  km,  and  0.093  m  at  90  km. 

At  altitudes  from  75  to  90  km,  the  space  vehicle  moves  in  the  flow  regime  close  to 
continuum  (the  Knudsen  number  changes  from  3.7  •  0~4  to  4.3  •  10”3).  The  drag  coefficient 
Op  at  altitudes  from  75  to  90  km  is  almost  independent  of  the  flight  altitude  and  the 
presence  of  chemical  reactions;  its  value  is  approximately  1.55.  The  heat-transfer  coeffi¬ 
cients  Ch  =  obtained  in  computations  with  and  without  chemical  reactions  are 

£o^  5* 

listed  in  Table  2.  It  is  seen  from  the  Table  that  the  heat-transfer  coefficient  decreases  by 
more  than  a  factor  of  3  at  an  altitude  of  75  km  and  by  a  factor  of  2  at  an  altitude  of  90 
km  if  chemical  reactions  are  taken  into  account. 

To  determine  the  efficiency  of  control  surfaces,  we  performed  three-dimensional  com¬ 
putations  of  spacecraft  aerodynamics  with  the  flap  deflected  by  0,  20°,  and  30°at  altitudes 
from  80  to  110  km  for  angles  of  attack  0  and  40°.  As  an  example,  Fig.  10  shows  the  pres¬ 
sure  held  and  the  surface  distribution  of  the  pressure  coefficient  at  an  altitude  of  80  km  at 
a  zero  angle  of  attack  and  40°.  It  is  seen  from  these  figures  that  the  control  hap  induces 
an  extremely  weak  perturbation  into  the  how  at  a  zero  angle  of  attack,  but  its  influence 
on  the  how  held  becomes  fairly  significant  as  the  angle  of  attack  is  increased  to  40°. 

The  DSMC  method  is  too  expensive  to  perform  a  multiparametric  study  of  aerother- 
modynamic  characteristics  of  the  space  vehicle  at  different  angles  of  attack,  different 
altitudes,  and  different  angles  of  hap  deflection. 

For  this  reason,  another  series  of  computations  was  performed  by  the  engineering  lo¬ 
cal  bridging  method.  To  estimate  their  accuracy,  the  results  of  these  computations  were 
compared  with  the  DSMC  results.  Figure  11  shows  the  aerodynamic  characteristics  as 
functions  of  the  angle  of  attack  for  a  hap  deflection  angle  of  30°.  The  curves  show  the 
results  for  the  angle  of  attack  equal  to  0  (1),  10  (2),  20  (3),  30  (4),  and  40°(5).  Dahsed 
curves  6  and  7  show  the  results  of  the  3D  DSMC  computations  for  the  angles  of  attack 
equal  to  40°and  0,  respectively.  There  is  no  curve  that  refers  to  the  zero  angle  of  attack 
in  the  Cn(H)  plot,  because  the  normal  force  is  close  to  zero  at  the  zero  angle  of  attack. 


Figure  9:  Overall  temperature  (left)  and  pressure  (right)  along  the  stagnation  line  at 
different  altitudes.  Axisymmetric  computation. 
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Altitude,  km  75  80  85  90 

Chemically  reacting  flow  1.30  •  10”2  2.11  •  10-2  3.38  •  10-2  6.18  •  10-2 
Chemically  inert  flow  4.30  •  10~2  6.22  •  10~2  9.47  •  10”2  1.36  •  lO^1 


Table  2:  Heat-transfer  coefficient  Ch  of  the  vehicle.  DSMC  results.  Non-deflected  flaps. 
Axisymmetric  model 


Cp  P,  Pa 

2 . 98e-03  1.48e+00 
1 . 19e-02  3 . 08e+00 
2 . 68e-02  4.82e+00 
4 . 77e-02  6.71e+00 
7 . 45e-02  8.79e+00 
1 . 07e-01  l.lle+01 
1 . 46 e-0 1  1 . 36 e+0 1 
1 . 91 e-0 1  1 . 64 e+0 1 
2 . 41e-01  1 . 96  e+O 1 
2 . 98  e-0 1  2 . 31 e+0 1 
3 . 61 e-0 1  2 . 71e+01 
4 . 29e-01  3 . 17e+01 
5 . 04 e-0 1  3 . 69e+01 
5 . 84e-01  4 . 31 e+0 1 
6 . 70 e-0 1  5 . 04 e+0 1 

7 .  63 e-0 1  5 . 91 e+0 1 

8 .  61 e-0 1  6 . 98 e+0 1 

9.  65e-01  8 . 31 e+0 1 
1 . 08e+00  1 . 00 e+0 2 
1 . 19e+00  1 . 23e+02 
1 . 31e+00  1 . 55e+02 
1 . 44e+00  2 . 03 e+0 2 
1 . 58e+00  2 . 83 e+0 2 
1 . 72e+00  4 . 43 e+0 2 


Cp  P,  Pa 

2 . 99e-03  1.51e+00 
1.20  e-0  2  3 . 16 e+0 0 
2 .  69e-02  4.94e+00 
4 . 79e-02  6.88e+00 
7 . 48 e-0 2  9.01e+00 
1 , 08  e-0 1  1 . 14 e+0 1 
1 . 47 e-0 1  1 . 39e+01 
1 . 91 e-0 1  1 . 68 e+0 1 
2 . 42e-01  2 . 00e+01 
2 . 99e-01  2 . 37 e+0 1 
3 . 62  e-0 1  2 . 78  e+0 1 
4 . 31 e-0 1  3 . 25 e+0 1 
06e-01  3 . 79e+01 
86 e-0 1  4 . 42 e+0 1 
73 e-0 1  5 .  l€e+01 
66e-01  6 .  06e+01 
65  e-0 1  7 . 15 e+0 1 
69e-01  8 . 52 e+0 1 
08 e+0 0  1 . 03 e+0 2 
20 e+0 0  1 . 2€e+02 
32 e+0 0  1 . 59e+02 
45 e+0 0  2 . 08 e+0  2 
58 e+0 0  2 . 90 e+0 2 
72 e+0 0  4 . 54 e+0 2 


Figure  10:  Pressure  held  and  surface  distribution  of  the  pressure  coefficient.  Altitude  80 
km.  Flap  deflection  angle  30°.  Angle  of  attack:  0  (left)  and  40° (right) 


The  moment  characteristic  differs  from  zero,  because  the  center  of  gravity  is  shifted  down 
on  0.156  m  from  the  spacecraft  axis.  The  curves  obtained  by  the  engineering  method 
describe  the  aerodynamic  coefficients  in  the  range  of  altitudes  between  120  and  60  km. 
The  difference  from  the  DSMC  results  is  approximately  5%  for  the  force  characteristics. 
The  accuracy  of  the  moment  characteristics  is  lower,  especially  at  an  angle  of  attack  of 
40°at  altitudes  of  100  and  110  km  (about  50%).  It  should  be  mentioned,  however,  that  the 
value  of  the  pitching  moment  coefficient  is  close  to  zero,  and  the  large  relative  difference 
actually  means  insignificant  absolute  difference  in  the  pitching  moment  coefficient.  The 
DSMC  computations  reveal  a  significant  effect  of  chemical  reactions  on  the  flow  fields 
around  the  spacecraft  and  on  the  heat-transfer  coefficient  at  altitudes  below  90  km.  Be¬ 
cause  of  considerable  changes  in  the  flow  structure  near  the  base  surface  of  the  spacecraft, 


Figure  11:  Axial  and  normal  force  coefficients  and  pitching  moment  coefficient  versus 
altitude  for  different  angles  of  attack.  Flap  deflection  angle  30°.  Angle  of  attack  0°(1), 
10°  (2),  20°  (3),  30° (4),  and  40° (5);  DSMC  results  for  the  angle  of  attack  equal  to  40° (6) 
and  0°(7). 
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there  is  a  several-fold  difference  in  the  values  of  the  pitching  moment  and  the  lift  force  of 
the  vehicle  with  deflected  control  surfaces,  which  were  obtained  with  and  without  chem¬ 
ical  reactions.  Simultaneously,  chemical  reactions  proceeding  in  the  gas  exert  practically 
no  effect  on  the  drag  coefficient.  The  changes  in  the  aerodynamic  characteristics  with 
variations  of  the  angle  of  attack  and  the  angle  of  flap  deflection  at  altitudes  from  80  to 
110  km  were  computed  by  the  DSMC  method.  The  aerodynamic  characteristics  of  the 
spacecraft  at  altitudes  from  120  to  60  km  were  also  calculated  by  the  engineering  local 
bridging  method.  Comparisons  with  the  DSMC  results  showed  that  the  axial  and  normal 
force  coefficients  are  calculated  rather  accurately  (within  5%).  The  results  obtained  can 
be  used  to  design  the  thermal  protection  system  of  the  space  vehicle  and  to  construct  its 
de-orbiting  trajectory. 

7  Prospects  for  the  DSMC  method 

In  conclusion,  let  us  mention  the  directions  for  the  DSMC  method  development  that 
we  believe  will  be  important  in  the  next  several  years. 

•  near-continuum  flows:  modeling  of  flows  at  relatively  high  Reynolds  numbers  (20,000 
and  higher)  using  both  conventional  DSMC  method  and  computationally  efficient 
hybrid  approaches 

•  real  gas  effects:  development  of  general  models  applicable  for  any  type  of  molecules 
that  would  accurately  predict  the  energy  transfer  processes  and  chemical  reactions 

•  numerical  accuracy  check:  use  of  available  and  development  of  new  criteria  that 
enable  one  to  evaluate  the  numerical  accuracy  of  results  obtained,  especially  in  the 
near-continuum  regime 

•  extensive  validation  of  models  and  algorithms:  comparison  of  DSMC  results  with 
experimental  data  and  contin-  uum  results  in  different  gas  dynamic  situations  for 
inert  and  reacting  gases 

•  efficient  parallel  algorithms  with  dynamic  load  balancing:  these  are  critical  for  mod¬ 
eling  of  computationally  intensive  two-  and  three-dimensional  problems 

8  Acknowledgments 


The  author  would  like  to  thank  the  research  team  of  the  Computational  Aerodynamics 
Laboratory  of  the  Institute  of  Theoretical  and  Applied  Mechanics,  Siberian  Branch  of 
the  Russian  Academy  of  Sciences.  Our  special  thanks  to  Yevgeniy  Bondar,  Alexander 
Kashkovsky,  Dmitry  Khotyanovsky,  Alexei  Kudryavtsev,  Gennady  Markelov,  Alexander 
Shevyrin,  Alina  Alexeenko  and  Pavel  Vashchenkov,  who  provided  results  of  computations 
and  participated  in  discussions  and  preparation  of  this  paper.  This  work  was  supported 
by  the  Russian  Foundation  for  Basic  Research  (Grant  No.  10-08-01203-a). 


18-20 


RTO-EN-AVT-194 


A  APPENDIX  1.  MAJORANT  COLLISION  FREQUENCY  SCHEMES 


A  Appendix  1.  Majorant  collision  frequency  schemes 


The  DSMC  method  is  traditionally  considered  as  a  method  of  statistical  simulation  of 
the  behavior  of  a  great  number  of  model  gas  molecules.  Usually  the  number  of  simulated 
particles  is  large  enough  (~  105  —  108),  but  this  is  extremely  small  in  comparison  with  the 
number  of  molecules  that  would  be  present  in  the  real  gas  flow.  Each  simulated  particle 
is  then  regarded  as  representing  an  appropriate  number  of  real  molecules. 

The  state  of  each  simulated  particle  is  characterized  by  its  coordinate  f  and  velocity  v. 
The  state  of  the  whole  system  of  N  particles  is  described  by  the  6 TV-dimensional  vector 
{R,  V}  =  {fi,  hi, ...,  Uv,  Av}.  The  evolution  of  such  a  system  can  be  represented  as  a 
jump-like  motion  of  a  point  in  the  6Y-dimensional  phase  space.  The  DSMC  method  can 
be  then  treated  as  statistical  simulation  of  the  6Y-dimensional  random  jump- like  process. 
For  such  a  simulation,  it  is  necessary  to  generate  the  trajectory  of  the  random  process. 
To  this  end,  it  is  needed  to  determine  a  technique  for  simulation  of  the  initial  trajectory 
point  of  the  random  process  (to,  Ro,  Vo),  transition  in  time  between  subsequent  collisions 
(free  motion)  from  the  state  (t0,  Ro,  Vo)  to  the  state  (U,  R i,  V0 ),  and  the  changes  in  particle 
velocities  after  collisions,  i.e.  the  transition  from  the  state  (U,  R\,  V0 )  to  (U,  Ri,  Vi).  Then 
transitions  from  state  1  to  state  2,  etc.  are  considered.  The  trajectory  is  terminated  when 
tn>  T  where  T  is  the  modeling  time. 

In  the  traditional  approach  of  constructing  numerical  schemes  of  the  DSMC  method 
[1],  the  description  of  procedures  for  trajectory  simulation  of  the  random  process  is  based 
on  physical  concepts  of  rarefied  gas  and  on  physical  assumptions  that  form  the  basis  for 
the  phenomenological  derivation  of  the  Boltzmann  equation. 

The  kinetic  theory  of  gases  makes  use  of  the  so-called  “master”  kinetic  equations  (MKE) 
[3,4]  which  describe  the  behavior  of  an  Y-particlc  gas  model  with  binary  collisions.  These 
linear  MKE  transform  to  the  nonlinear  Boltzmann  equation  when  N  — *  oo  and  molecular 
chaos  conditions  are  satisfied  (see,  e.g.,  [5]  ).  Since  a  finite  number  of  simulated  particles 
is  used  in  numerical  simulation,  it  is  natural  to  utilize  directly  these  MKE  equations  for 
constructing  numerical  schemes  of  the  DSMC  method.  Such  an  approach  to  the  derivation 
of  numerical  schemes  of  the  DSMC  method  was  presented  in  [2,33,37]. 

First,  let  us  consider  a  spatially-uniform  rarefied  gas  flow,  and  present  the  main  steps 
of  constructing  numerical  schemes  directly  from  the  MKE.  In  this  case,  it  is  possible  to 
use  the  Kac  master  equation  [3]: 

i)t  J'xU-  V)  =  /  d6ii  /  MMu  -  Vj\{fN(t,  V'j)  -  fN(t ,  V )},  (1) 

where  V  =  (hi, . . .  ,vn)  is  a  3Y-dimensional  vector  of  particle  velocities;  Vt(  =  (hi, . . . , 
v(, . . .  ,Vj, . . . ,  hjv);  bij  and  etJ  are  the  impact  parameters;  (h',h')  and  (vi,Vj)  are  the  pre- 
and  post-collisional  velocities  of  i ,  j  molecules,  n  is  the  number  density,  and  the  sum  over 
i  <  j  means  the  summation  over  N(N  —  l)/2  collision  pairs. 

This  is  a  linear  integro-differential  equation  that  describes  the  time  behavior  of  the 
Y-particle  distribution  function  fN(t,V),  f  fjy(t,V)dV  =  1. 

In  accordance  with  the  general  theory  of  Monte  Carlo  methods  [39,40],  the  following 
approach  is  used:  the  transition  from  an  integro-differential  form  of  the  master  kinetic 
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equation  to  a  linear  integral  equation  whose  probability  treatment  serves  as  the  basis  for 
constructing  an  appropriate  random  process  of  direct  simulation. 

Using  a  function  w  -  the  probability  density  of  transition  of  a  pair  of  particles  from 
(u',h')  to  (vi,Vj),  eqn  (1)  is  written  as 

r\  s* 

V)  +  u(V)fN(t,  V)  =  ^  Vij)w(yv'i,  v'j  ->■  v^dv-dv'j,  (2) 

i<j  d 


where  u(V)  is  the  total  collision  frequency 


N 


l4V)  =  —  Y  I  wWiv'j  ->■  vh  Vj)dVidvj  =  —  Y  I  a  ( 9ij ,  £ij)  h  (vi  +  v3  -v[-  v'3)  x 


i<j 


1=3 


xJi 


(Vi  -  Vj)2  -  (v'i  -  v'j)- 


dv'idv'j  =  4  Y  a49ij)9ij  < 


oo 


(3) 


i<j 


and  g  is  the  relative  collision  velocity,  (Jt{9ij)  and  (j(gt-j,  QQ  are  the  total  and  differential 
collision  cross-sections,  and  Xij  is  the  deflection  angle. 

The  total  collision  frequency  (3)  determines  the  time  of  the  next  collision  in  the  system 
and  depends  on  velocities  of  all  particles,  and  it  is  necessary  to  recalculate  its  value  after 
each  collision.  This  is  rather  time-consuming  if  the  number  of  particles  N  is  large,  since 
the  summation  is  performed  over  all  possible  N(N  —  l)/2  collision  pairs. 

Let  us  introduce  the  majorant  collision  frequency  [33] 


N(N  —  1)  r  , 

On  =  — - [9&t\9 Jjmax  >  V{V)- 


(4) 


Then  we  add  to  the  both  sides  of  eqn  (2)  the  corresponding  sides  of  the  equality 

K  -  v)  =  ^  f  M*’ i ~  AMs't,)}  W  -  v')dv' 


i<j 


and  join  the  right-side  integrals  to  have 


d 


n 

N 


-fN(t,V)  +  vmfN(t,V)  = 

Y  I  [[9Vt(g)]max  -  g'ijVtig'ij)]  SiVi-v'JS^Vj-Vjj+W^Vj  ->■  Vi}Vj)  j>dh'du'. 

(5) 


i<j 


Here  <5  denotes  Dirac’s  delta-function. 

Let  us  introduce  the  function  V)  as 


fN(t,  V)u{V)  =  /  K2(t'  t\V)iP(t',  V)dt', 


where  K2  is  the  kernel 


K-2(t'  — >  t\V)  =  9{t  —  t')u(V)  exp  {— v(V)(t  —  t')} 
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and  9  is  the  Heaviside  function,  then  eqn  (5)  may  be  transformed  to  an  integral  form 


K^{t',V')dV'dt'  +  5(t)fUV) 


with  the  kernel 

Here 


K2 1  =  K2 (P  tlV^K^V'  ->■  V). 

K2(t'  ->■  t\V')  =  9{t  -  t')vm  exp{— z/TO(f  -  P)}, 

2 


Kl(V'^V)=Y. 


i<j 


N(N-l) 


x 


x 


[gvt{g)]mc 


5(vi-v'i)8(vj-v'j)  + 


g'ijMg'ij)  wWvVj  ->•  Vi,Vj 


\g°t(g)\ , 


g'iMg'ij] 


(6) 


m=l 


The  probability  treatment  of  eqn  (6)  enables  one  to  construct  a  random  process  that 
describes  the  behavior  of  the  Y-particlc  gas  model.  To  calculate  the  initial  trajectory 
point  of  this  process,  the  free  term  from  eqn  (6)  is  used  as  the  probability  density,  while 
the  kernel  K2i  yields  the  probability  density  of  the  transition  from 


the  state  (P,  V')  to  the  state  (t,V). 


This  transition  is  modeled  sequentially  from 

the  state  (P,  V')  to  (t,  V') 
in  accordance  with  the  distribution  density  A'2(P  — »  t\V'). 

Then,  in  compliance  with  the  kernel  K \  ( V'  — >  V)  the  transition  (t,  V')  — >  (t,  V)  occurs. 

For  this  purpose,  a  collisional  pair  (i,j)  is  uniformly  chosen  from  N(N  —  l)/2  pairs 
(this  is  evidenced  by  the  presence  of  the  factor  2/N(N  —  1)  in  the  kernel  Ad).  A  collision 
of  this  pair  occurs  with  the  probability 

p_  g'Ma'ij) 

[gvtig)]  max 

and  with  the  probability  (1  —  P)  the  velocities  will  not  change  (this  is  evidenced  by  the 
presence  of  the  product  of  two  delta-functions  in  the  first  term  of  the  kernel  Ah),  i.e.  a 
fictitious  collision  occurs.  The  product  of  delta-functions  in  K\  (the  last  factor)  shows  that 
after  a  collision  of  the  pair  (i,j)  velocities  of  other  particles  remain  unchanged,  whereas 
velocities  (h',h')  for  real  collisions  are  replaced  by  post-collisional  velocities  ( Vi,Vj ). 
Collision  impact  parameters  are  selected  with  the  probability 

W(v'i:v'j  -»  Vi,Vj) 

g’ijVtWij ) 

and  the  presence  of  two  delta-functions  S3  and  <5i  in  w  (see  (3))  shows  that  new  velocities 
are  calculated  in  accordance  with  momentum  and  energy  conservation  laws. 
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Thus,  all  required  procedures  for  trajectory  simulation  of  the  random  process  are 
specified,  and  the  numerical  scheme  of  the  DSMC  method  is  constructed  for  a  spatially 
uniform  case.  Its  computer  cost  is  proportional  to  the  number  of  simulated  particles  [2], 
Usually  in  problems  of  rarefied  gas  dynamics  it  is  required  to  determine  gas  parameters 
at  the  time  moment  A,  which  have  the  form 

4(4)  =  I  h(v)fi(tkv)dv  =  I  H(y)fNdV , 

where  h(v)  is  the  function  of  velocity,  and 

To  calculate  functionals,  such  as  4(4),  it  is  possible  to  use  the  estimates  known 
from  the  general  theory  of  Monte-Carlo  methods  [40].  In  particular,  a  counterpart  of  the 
non-biased  absorption  estimate  has  the  form 

£,/,  =  H(VS),  s  =  max{n  :  tn  <  L,  n  =  0, 1, ...}. 

Approximate  values  Ih(tk )  of  the  functionals  4(4 )  are  given  by 

L 

ih(tk)  =  L~1J2U  0, 

i=i 

where  L  is  the  number  of  independent  Y-particle  trajectories  of  a  random  process.  Fol¬ 
lowing  [33],  it  can  be  shown  that  the  mathematical  expectation  is 

£'[&/’]  =  4(4), 

and  the  variance  of  random  variable  is 

Var  [At,]  =  Y_1  j  j"  h 2  (u)  /j  (t,  v)  dv  -  l‘l  (4)  j  + 

H — j —  j  h  (ui)  h  (v2)  {f'2  ( t ,  vi,  v2)  -  fi  (t,  vi)  f  (t,  v2 )}  dvidv2.  (7) 

It  is  known  [39]  that  the  variance  of  absorption  estimate  is  limited  and,  by  the  central 
limit  theorem,  the  following  inequality  is  valid: 

L 

<3 

i=i 

It  is  commonly  believed  that  the  statistical  error  of  the  DSMC  method  is  determined 
only  by  the  overall  volume  of  the  sampling,  (N  ■  L).  It  follows  from  (7)  that  under  the 
condition  of  molecular  chaos  (f2  =  fifi)  the  second  term  vanishes,  and  the  statistical 
error  is  proportional  to 
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For  finite  N,  however,  there  is  always  a  statistical  dependence  of  particles  (see  [33]  for 
more  detail);  and,  hence,  the  statistical  error  is  proportional  to 

i/Vl. 


To  analyze  the  relationship  between  simulation  results  and  the  solution  of  the  Boltz¬ 
mann  equation,  it  is  necessary  to  consider  the  kinetic  equation  for  the  one-particle  distri¬ 
bution  function 

A  (Mi)  =  J  fN{t,VN)dv2...dvN. 

If  Kac  MKE  equation  (1)  is  integrated  over  (N  —  1)  variables  v2) ...,  vn,  then  we  obtain 
a  kinetic  equation  for  the  one-particle  distribution  function 


d_ 

Ft 


A 


N  —  1 


N 


n 


bdbx 


x  J  dv2 \vi  -  u2|(/iK)/i(u2)  -  fi(vi)fi(v2))+ 


N  —  1 


N 


n 


dv2\vi  -  v2\(g'2  -  g2), 


(8) 


where  g2  —  f2  —  fifi-  In  the  general  case  the  solution  of  equation  (8)  differs  consider¬ 
ably  from  the  solution  of  the  Boltzmann  equation,  since  it  depends  on  the  variation  of 
the  two-particle  correlation  function  g2(t,v i,v2).  Examples  of  the  influence  of  statistical 
dependence  of  particles  on  the  simulation  results  were  presented  in  [33].  Equation  (8)  is 
transformed  into  the  Boltzmann  equation  [41]  when  N  — >  oo  and  the  molecular  chaos 
hypothesis  is  valid.  Strictly  speaking,  only  under  these  conditions  the  results  of  statistical 
simulation  precisely  correspond  to  the  solution  of  the  Boltzmann  equation. 

It  is  shown  therefore,  that  the  use  of  the  MKE  allows  one  to  obtain  not  only  a  prob¬ 
abilistic  procedure  for  simulation  of  the  random  process  and  statistical  estimates  for  cal¬ 
culating  the  gas  parameters,  but  also  to  assess  the  statistical  error  of  results.  In  addition, 
the  relationship  between  the  results  of  simulation  of  the  Y-particlc  gas  model  and  the 
solution  of  the  Boltzmann  equation  has  become  more  transparent. 

The  next  step  is  to  illustrate  how  the  approach  described  can  be  extended  to  a  spa¬ 
tially  nonuniform  case.  The  traditional  approach  [1]  employs  the  discretization  of  spatial- 
temporal  evolution  of  the  Y-particle  gas  model.  A  continuous  process  of  motion  and 
collisions  of  molecules  is  uncoupled,  the  computational  domain  is  divided  into  cells,  and 
the  following  processes  are  sequentially  simulated  at  each  step  At: 


•  spatially  uniform  relaxation  of  the  gas  in  each  cell; 

•  free-molecular  movement  of  simulated  particles  over  the  distances  appropriate  to  At 
with  regard  for  boundary  conditions. 

A  different  approach  was  proposed  by  Ivanov  and  Rogasinsky  [2]  who  constructed 
numerical  schemes  for  simulating  a  6Y-particlc,  continuous  in  time,  random  process  of 
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spatially  nonuniform  evolution  of  a  system  of  N  particles.  These  schemes  were  derived 
directly  from  the  Leontovich  master  kinetic  equation  [4]: 


N 


dfN  +  ^  VidI^-  fN  =  S(n  -  rj)  I  { f'N  -  fN}\vi  -  Vjlbijdbijdeij.  (9) 


dt 


1=1 


i<j 


Here  J'm  =  fn(t,  R,  V),  f'N  —  fn(t ,  R,  U/)  are  the  Y-particle  distribution  functions,  and 


fNdVdR  =  1. 


The  relationship  between  this  equation  and  spatially  nonuniform  Boltzmann  equation 
was  studied  in  [5]. 

The  presence  of  the  delta-function  in  collision  integral  (9)  indicates  the  collision  “lo¬ 
cality”,  i.e.  identical  coordinates  of  colliding  particles,  like  in  the  Boltzmann  equation. 
For  numerical  simulation  one  has  to  perform  the  regularization  of  collision  integral  (9), 
i.e.  introduce  “smeared”  particle  interactions.  Using  the  function  wp,  the  collision  integral 
in  eqn  (9)  takes  the  form 


Jn  — 


N  r 

J  ( fN-fN 

l=J 


wpdv{dvj, 


where 

is  the  probability  density  of  transition  of  the  pair  (i,j)  from  the  state  (v'^Vj)  to  the  state 
(vi,  Vj)  with  fixed  values  of  particle  coordinates  (fl,  f3 )  and  regularization  parameter  p.  In 
this  case, 


wp  (v[,  Vj  ->■  Vi,  Vj\ri,  rj,  p j  ->•  5  (v i  -  rj)  w  (v^v'j  ->■  v^  v^j  when  p  -»  0 
Let  us  consider  the  following  regularization 


w, 


(yi,Vj  ->•  v^Vjlr^rj,^  =  h(ri,rj)w  (v'^v'j  ->  v^v^j  , 


h(ri,  rj) 


w0\  if  \ri~rj\  <  p, 
0,  if\fi-fj\>p, 


=  g  71 P3- 

The  parameter  p  defines  the  dimension  of  the  “interaction  region”  of  particles,  and  its 
value  depends  on  the  required  accuracy  of  calculation  of  flow  parameters. 

Then,  the  total  collision  frequency  of  an  Y-particle  gas  model  is  as  follows 


N 

isp(R,  V)  =  ^  h(r i,  fj)gijat(gij)  <  oo,  (10) 

i<j 
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and  it  depends  on  the  particle  coordinates  and  velocities.  As  in  the  case  of  a  uniform 
gas  model,  the  calculation  of  this  collision  frequency  requires  the  summation  of  N(N  — 
l)/2  pairs  with  checking  the  distance  between  the  particles,  which  is  time-consuming. 
Therefore,  here  it  is  again  profitable  to  majorize  the  collision  frequency,  and  contrary  to 
(4)  not  only  with  respect  to  the  particle  velocities,  but  also  to  the  coordinates  of  colliding 
particles. 

Let  us  choose  the  majorant  collision  frequency  in  the  form 

N(N-l)  i  -  - 

Dn  =  - 2 - ~U0  [gOt{g)]max  >  »{R,  V ),  (11) 

and  transform  the  Leontovich  equation,  similar  to  the  uniform  case,  to  the  following  form: 


dfN 

dt 


dfN 


N 

+  +  VmfN  = 

OT  i 

i=l 


=  fN{h(.ri,rj)uJ01w(v,i,v'j  ^Vi,Vj)+ 

i<j  J 

+  ({gvtig^maxUo1  ~  h[f  i,  r^u^gijatigij))  x  5{v[  -  u')<5(b'  -  v')}du'c%.  (12) 

The  collision  integral  here  consists  of  two  parts  corresponding  to  real  (first  term)  and 
fictitious  (second  term)  collisions. 

The  use  of  such  a  majorant  collision  frequency  um  allows  one,  therefore,  to  trans¬ 
form  the  nonuniform  problem  under  consideration  to  a  spatially  uniform  problem  with  a 
constant  collision  frequency. 

If  eqn  (12)  with  appropriate  initial  and  boundary  conditions  is  transformed  to  a  linear 
integral  equations  (see  details  in  [2]),  then  the  probability  treatment  of  the  kernel  and  free 
term  of  this  equation  allows  for  formulating  a  new  scheme  of  direct  statistical  simulation 
of  spatially  nonuniform  rarefied  gas  flow  with  continuous  time. 

Below,  the  salient  points  in  simulation  of  a  6Y-dimensional  trajectory  of  the  random 
process  of  the  transition  of  an  Y-particle  system  from  the  state  (P,  R',  V')  to  (t,  R,  V)  are 
briefly  described. 

We  shall  not  dwell  here  upon  the  simulation  of  particles  entering  the  computational 
domain  and  the  interaction  of  particles  with  the  body  surface.  All  details  can  be  found 
in  [2], 

The  time  of  the  next  collision  t  has  the  probability  density  distribution 

z/mexp  -  t')}  , 

while  the  probability  density  of  the  transition  of  a  system  of  N  particles  from  the  state 
(P,  R' ,  V')  to  (p  R,  V')  is  the  delta-function 

5{R-  R1  -V\t-t')). 

This  means  that  all  simulated  particles  move  during  the  time  t  —  P  from  point  R'  to  R 
with  velocities  V' . 
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An  examination  whether  the  collision  is  real  is  performed  as  follows: 

if  |  f,  —fj  |  >  p  then  the  collision  is  fictitious,  else  the  collision  is  real  with  the  probability 

Hij^t !//'.; ) 

max 

and  fictitious  with  the  complementary  probability. 

Thus,  all  probabilistic  procedures  necessary  for  the  numerical  generation  of  the  tra¬ 
jectories  of  a  6Y-dimensional  random  process  of  evolution  of  a  spatially  nonuniform  N- 
particle  gas  model  are  described. 

For  regularization  of  the  collision  integral  in  eqn  (9)  it  is  also  possible  to  use  a  usual 
approach  of  the  DSMC  method:  the  computational  domain  is  divided  into  non-intersecting 
cells  dk,  so  that 

M 

^  dk  =  V0, 

k= 1 

and  a  collision  may  occur  only  between  particles  that  belong  to  the  same  cell.  Then  the 
majorant  collision  frequency  has  the  form 


m 


N(N  —  1) 


[gvtia)} 


maxdmin ,  dmin  nim  (dk) 


(13) 


A  collision  of  the  pair  {i,  j)  randomly  chosen  from  all  N(N  —  l)/2  pairs  is  fictitious  if  the 
particles  i  and  j  are  in  different  cells. 

The  use  of  the  majorant  collision  principle  for  spatially  nonuniform  case  with  the 
regularization  of  the  type  (11)  or  (13)  allows  one,  therefore,  to  obtain  a  new  exact  time- 
continuous  scheme  of  the  DSMC  method. 

A  salient  point  that  differs  this  scheme  from  the  conventional  DSMC  method  is  the 
transfer  of  all  particles  after  each  collision.  Certainly,  this  is  extremely  expensive;  there¬ 
fore,  only  examples  of  using  this  algorithm  for  ID  rarefied  gas  dynamic  problems  (shock 
wave  structure  and  Couette  flow)  are  now  available.  However,  the  mere  fact  that  it  is 
possible  to  simulate  rarefied  gas  flows,  avoiding  uncoupling  of  a  continuous  process  of 
molecular  motion  and  collisions,  is  of  principal  importance.  Possibly,  this  approach  can 
be  used  in  future  when  the  computer  performance  significantly  increases. 

It  is  natural  to  consider  an  approximate  simulation  scheme  with  separated  and  se¬ 
quentially  calculated  collisions  and  motions  of  particles. 

Let  us  introduce  discrete  time  instants  tn  =  n  ■  At  and  choose  a  step  At  such  that  for 
each  rm,  provided  )T)m  rm  <  At,  the  following  condition  is  valid: 


vp  ( R  +  Tm  ■  V,  V\p)  «  up  (R,  V\p) 


This  means  that  the  change  in  collision  probability  of  pairs  due  to  particle  displacement 
is  ignored  during  the  time  At.  Then  an  approximate  simulation  scheme  is  as  follows:  the 
time  between  consecutive  collisions  is  chosen  on  the  basis  of  the  probability  density 


z/mexp{-z/m(r)} . 


K 

J2  T m  < 
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then  a  collision  (real  or  fictitious)  occurs.  All  particles  are  displaced  at  a  distance  pro¬ 
portional  to  At  one  time,  at  the  end  of  the  step  At. 

Numerical  schemes  with  a  continuous  and  discrete  time  that  make  use  of  the  majorant 
collision  frequency  with  the  regularization  of  types  (11)  or  (13)  do  not  employ  the  sorting 
of  particles  over  cells  and,  therefore,  can  be  called  free  cell  schemes.  The  determination 
of  distributed  and  total  aerodynamic  characteristics  can  be  conducted  directly  using  only 
the  collisions  of  a  particle  with  the  body  surface.  The  calculation  of  flowfields  requires  a 
grid  for  sampling  the  gas  parameters,  but  this  grid  is  not  related  to  collision  simulation. 

The  value  of  the  free  cell  majorant  can  be  reduced  by  sorting  the  particles  over  the 
ppIIg  Thpn 

r  fv  Nk{Nk  -  1)  \gat{g)}kmax 
Vrn~2^  vm-2^  2  4 

k= 1  k= 1 

where  is  the  number  of  molecules  in  the  /c-th  cell.  Obviously, 

, , cell  ^  , ,  freecell 

and,  hence,  the  computational  cost  of  the  cell  scheme  is  smaller  than  that  of  the  free  cell 
scheme.  In  the  cell  scheme,  however,  it  is  necessary  to  take  into  account  the  cost  of  particle 
sorting  over  cells.  This  sorting  is  very  simple  in  the  case  of  a  rectangular  Cartesian  grid, 
but  with  using  more  complex  adaptive  body-fitted  grids  its  cost  drastically  increases. 

When  using  the  cell  regularization  with  the  majorant  (14),  having  determined  the  time 
of  the  next  collision  Tlm  at  the  step  At,  one  has  to  determine  the  number  of  the  cell  in 
which  this  collision  occurs.  The  cell  number,  k,  is  selected  with  the  probability  v^n/vm. 
Then  a  pair  (i,j)  is  uniformly  chosen  in  the  k-th  cell  from  particles.  The  collision  is 
real  with  the  probability 

dij&tjgij) 

[g°t{g)}  max 

and  fictitious  with  the  complementary  probability. 

The  mean  number  of  collisions  in  a  system  of  N  particles  during  the  time  At  is  vmA t, 
and  each  cell  is  chosen  with  the  probability  zA/ vm.  Consequently,  during  the  time  At 
each  cell  is  chosen  z/^A t  times  in  the  average.  Then  for  z/^A t  3>  1  there  is  no  need  to 
sample  the  cell  number,  and  the  usual  sequential  computation  of  collisions  over  all  cells 
can  be  employed. 

The  scheme  that  uses  the  cell  regularization  with  particle  sorting  and  probabilistic 
choice  of  the  cell  number  where  a  collision  will  occur  can  be  called  a  cell  majorant  collision 
frequency  scheme. 

Our  experience  of  using  the  cell  and  free  cell  schemes,  which  has  been  accumulated 
up  to  now,  shows  that  these  schemes  should  not  be  set  in  opposition.  Vice  versa,  it  is 
desirable  to  use  them  together.  For  this  purpose,  the  computational  domain  is  divided 
into  uniform  Cartesian  cells  with  the  linear  dimension  of  the  order  of  Aqo,  and  the  cell 
scheme  is  used  there.  In  flow  areas  where  the  local  mean  free  path  strongly  decreases  and 
in  the  vicinity  of  the  body  surface,  the  free  cell  scheme  is  used  with  a  small  regularization 
parameter  in  order  to  increase  the  spatial  collision  resolution. 

Now  a  few  words  about  the  influence  of  the  finite  number  of  simulated  particles  on 
the  results  of  computations.  In  the  nonuniform  case,  the  main  criterion  that  allows  for 
a  practical  verification  of  the  presence  of  a  statistical  dependence  between  the  simulated 


RTO-EN-AVT-194 


18-29 


A  APPENDIX  1.  MAJORANT  COLLISION  FREQUENCY  SCHEMES 


particles  is  the  relative  number  of  repeated  collisions.  Here,  repeated  collisions  are  colli¬ 
sions  between  the  same  pair  of  particles  during  their  lifetime  in  the  computational  domain. 
This  number  is  directly  related  to  the  number  of  particles  N\  in  a  volume  with  linear  size 
of  the  local  mean  free  path.  If  N\  >  1,  we  can  definitely  say  that  the  simulation  results 
are  close  to  the  solution  of  the  Boltzmann  equation. 

This  is  certainly  only  an  intuitive,  experience-based  criterion.  A  more  detailed  study 
of  the  influence  of  statistical  dependence  with  the  use  of  autocorrelation  functions  and 
the  number  of  repeated  collisions  can  be  found  in  [42], 

The  extension  of  the  majorant  collision  frequency  schemes,  directly  obtained  from  the 
MKE,  to  the  case  of  multi-species  chemically  reacting  mixtures  is  straightforward.  It 
implies  the  change  in  cross-sections  of  the  corresponding  inelastic  processes  and,  hence, 
the  change  in  the  collisional  algorithm  only  in  its  part  that  refers  to  collision  mechanics. 
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A  Appendix  2.  Gas/surface  interaction  models 


Calculation  of  aerodynamic  forces  by  the  DSMC  method  requires  setting  the  veloc¬ 
ity  distribution  functions  for  molecules  reflected  from  the  object  surface.  A  simplified 
representation  of  this  function  is  used,  which  must  take  into  account  the  major  features 
of  gas-surface  interaction  which  have  been  revealed  in  experiments  (force  action,  angu¬ 
lar  distribution  of  molecules  escaped),  on  the  one  side,  and  it  must  be  simple  enough  in 
application,  on  the  other  side. 

The  main  requirement  of  the  surface  properties  definition  is  to  introduce  for  every 
element  of  the  space  object  the  following  surface  properties: 

1.  Wall  temperature 

2.  Scheme  of  gas-surface  interaction  (Maxwell  or  Nocilla)  model,  accommodation  and 
interaction  coefficients. 

If  a  space  object  consists  of  several  hundred  or  thousand  elements  it  is  not  convenient 
for  user  to  introduce  or  edit  the  properties  of  each  element.  So  a  procedure  was  found  to 
give  the  user  a  fast  possibility  for  introducing  and  editing  with  a  maximum  flexibility: 

1.  Setting  surface  properties  for  the  geometric  model.  For  each  compound  and  element 
a  list  of  surface  properties  can  be  assigned. 

2.  Introducing  and  checking  the  surface  properties  of  a  compound.  The  system  shows 
the  actual  content  of  this  compound  and  to  each  element  surface  properties  can  be  as¬ 
signed. 

3.  To  check  the  correct  settings  of  the  surface  properties  the  system  provides  a  visu¬ 
alization  tool. 

One  of  the  oldest  and  most  widely  used  distribution  functions  for  the  reflected  molecules 
is  due  to  Maxwell  [43]. 

It  is  constructed  on  the  assumption  that  a  fraction  1  —  a  of  the  molecules  is  reflected 
from  the  surface  in  a  specular  fashion,  while  the  fraction  a  is  re-emitted  diffusely,  with 
the  Maxwellian  distribution: 


fr(r,£r) 


(1  -  a)  fi  (r,  £r  -  2(£r  •  n)n)  +  a  vr  3/2  c'r  3  exp 


where  is  the  reflected  velocity  ft  is  the  external  normal  of  the  surface  in  point  r,  c'r  is 
the  most  probable  speed  at  the  temperature  Tr. 

The  parameter  a  is  called  the  tangential  momentum  accommodation  coefficient.  Ev¬ 
idently,  for  completely  specular  reflection,  o  =  0,  while  for  complete  diffuse  reflection, 
cr  =  l. 

This  implementation  of  the  specular- diffuse  scheme  implies  either  specular  or  diffuse 
reflection  of  particles.  Therefore,  at  first,  each  particle  is  sampled  to  experience  a  par¬ 
ticular  type  of  reflection.  If  the  pseudo-random  number  is  smaller  than  <7,  the  particle 
experiences  diffuse  reflection;  otherwise,  it  experiences  specular  reflection. 

In  the  case  of  specular  reflection,  the  velocity  component  normal  to  the  surface  changes 
its  sign,  and  the  tangential  components  of  velocity  remain  unchanged. 
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The  velocity  of  diffusely  reflected  molecules  is  sampled  in  the  local  polar  coordinate 
system,  directed  along  the  outer  normal  of  the  to  surface  of  the  vehicle,  in  the  following 
manner: 


I  &■  |=  4  a/-  ln(i?iR2); 
cos(0)  =  v^s; 

0  =  2  7ri?4, 

where  |  £r  |  -  module  of  the  velocity  of  reflected  particle;  0,0-  polar  and  azimuthal  angles 
and  R  \ .  R> .  R:> .  R4  are  random  numbers; 

Following  ref.  [43]  an  energy  accommodation  coefficient  oe  =  Ife— is  used.  When 
the  accommodation  coefficient  Oe  =  0  then  the  molecules  don’t  exchange  energy  with  the 
wall,  and  Oe  =  1  if  the  incident  molecules  reached  thermal  equilibrium  with  the  wall. 

If  we  write  the  energy  of  a  molecule  in  terms  of  its  velocity,  Eq.  A  takes  the  form: 

42-K2, 

E  Q  -  c'r2  ’ 

accordingly,  the  value  of  velocity  with  which  the  particles  are  reflected  is  written  as 

K  =  y/(l  —  Te)42  +  °ec'1 

Thus,  for  each  diffusely  reflected  molecule,  the  absolute  value  of  its  sampled  velocity 
|  |  should  be  corrected  by  the  coefficient 

k  =  lr  =  \Z(!  —  °E)tf/c'r  +  (Te, 

and  hence 


|  |=  kc'r  \J  —  \n(RiR2). 

The  distribution  function  [44-45]  for  the  velocities  of  reflected  particles  which  is  often 
used  in  practice  has  the  following  form: 


fr  Tl  r  ( 


1  m 


7 r  2kTr 


exp 


m 


2  kTr 


Z-ur 


1 2 


1  m  \  2 

=  n'Rl Wr)  ^ 


m 

2kTr 


t-Sr 


=  nr\- 


1  m 


ix  2kTr 


exp<  — 


m 

2kTr 


ty  ~  Sn 


m 

2kXr 


L-sr 


m 

2kTr 


G  'rr,z 


Here  m  and  $,  are  the  mass  and  velocity  of  molecules  reflected  from  surface,  and 
parameters  of  particle  number  density  nr ,  temperature  Tr  and  speed  Ur  depend  on  the 
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velocity  vector  of  a  molecule  incident  upon  the  surface,  as  well  as  on  its  energy,  surface 
characteristics  and  total  mass  flow  of  incident  molecules.  These  dependences  transform 
to  those  of  five  scalar  parameters  of  function  fr  on  the  base  of  available  experimental  data 
and  mass  conservation  law. 

The  experiments  [44]  showed  that  the  vector  of  the  mean  reflected  velocity  is  located  in 
the  same  plane  as  the  vector  of  the  incident  velocity.  Thus,  if  we  construct  the  right-hand 
Cartesian  coordinate  system  at  the  point  of  particle  reflection  so  that  the  x  axis  is  the 
projection  of  the  incident  velocity  vector  onto  a  plane  tangential  to  the  surface,  the  y  axis 
coincides  with  the  external  normal  to  the  surface,  then  the  side  component  of  reflected 
velocity  is  STr-tZ  =  0. 

The  four  parameters  of  this  function  Snr,  STt .)X,  Tr,  nr  are  found[44,45]  as  follows: 

29 

Snr  =  0.1  —  0.65 - , 

71 


q  _  7~r,x  MSnr) 

^ rr,x  '  j  /  o  \i 

Pr  Jl\&nr) 


k 

m 


2 —Tr  = 


f— V 

rJl(Snr)] 

\  Qmr  / 

Qmr  (  2  Tr  )  d  ]  ( Snr  ]) , 
m 


T, 


r,  x 


40  —  7r\ 

2vr  J 


£  sin{9)  qmi 


Pr  = 


Ti  kTw 
2  m 


Qmi 


Jiit)  =  -=  qexp  [-(q  -  t)2]  dq  =  — =  e  +  ^Txt  (1  +  erf(t)) 

"  0  V27T  L 


1  1 

Mt)  =  —j=  /  d2  exp  [~(q  -  t )2]  dq  = 


7T 


L 

erf[t)  =  —=  [  e~s2  ds, 


te~t2  +  \pn  (-  +  t2)  (1  +  erf(t)) 


TV 


Qmr  =  Qmi  =  mq^  qi  =  1  [m  s  . 


Here  9  is  the  angle  between  the  incident  velocity  vector  and  internal  normal  to  the  surface, 
Tw  is  the  surface  temperature,  qmr  and  qmi  are  the  fluxes  of  reflected  and  incident  molecules, 
which  is  used  to  normalize  the  distribution  function  fr  in  the  Monte-Carlo  method. 

The  values  of  model  parameters  an,  bn,  aT,  and  bT  depend  on  the  type  of  material,  and 
the  following  values  were  obtained  [44,45]  based  on  experimental  data: 
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Materials 

&n 

|  K 

aT 

bT 

Ceramics 

0.06 

0.05 

-0.03 

0.26 

White  enamel 

0.07 

0.10 

-0.01 

0.16 

Black  enamel 

0.065 

0.08 

-0.03 

0.11 

Scree-vacuum  glass  Cloth 

0.07 

0.03 

-0.05 

0.11 

Anodized  aluminium- magnesium  Alloy  (AMA) 

0.10 

0.10 

-0.06 

0.04 

Chemically  polished  AMA 

0.11 

0.18 

-0.19 

-0.04 

Glasses 

0.17 

0.09 

-0.10 

-0.08 

Enamels 

0.07 

0.06 

-0.01 

0.12 

Sampling  of  velocity  with  which  a  particular  particle  is  reflected  involves  the  proba¬ 
bility  density  of  reflected  molecules,  which  is  obtained  from  the  distribution  function  and 
is  written  in  the  following  form: 


p  — 

■L  r.T. 


\fTX\J% 


exp< 


m 

2kTr 


€x  -  ST 


F  = 

1  r,y 


V^VTrMSnr)  \2kT, 


m 


Cy  exp<  - 


m 

2  kTr 


k,y  Snr 


F  — 

1  r,z 


\pK\J% 


exp< 


l  n  2  > 

m  \ 2 
2kT,  1  * z 


The  velocity  £,y  is  simulated  by  the  Neumann  method,  and  the  velocities  £x  and  are 
simulated  on  the  basis  of  the  relations 


=  \J -  ln(f?i)  cos(2n  R2)  +  Srr, 


x 


-j  =  yj—  ln(.Ri)  sin(2ir  R2) 

C'r 

where  R\  and  R2  are  pseudo-random  numbers. 
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